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Abstract 

The recent two-dimensional (2-D) parametric formulation of the high fidelity generalized method of 
cells (HFGMC) reported by the authors is generalized for the micromechanical analysis of 
three-dimensional (3-D) multiphase composites with periodic microstructure. Arbitrary hexahedral subcell 
geometry is developed to discretize a triply periodic repeating unit-cell (RUC). Linear 
parametric-geometric mapping is employed to transform the arbitrary hexahedral subcell shapes from the 
physical space to an auxiliary orthogonal shape, where a complete quadratic displacement expansion is 
performed. Previously in the 2-D case, additional three equations are needed in the form of average 
moments of equilibrium as a result of the inclusion of the bilinear terms. However, the present 3-D 
parametric HFGMC formulation eliminates the need for such additional equations. This is achieved by 
expressing the coefficients of the full quadratic polynomial expansion of the subcell in terms of the side or 
face average-displacement vectors. The 2-D parametric and orthogonal HFGMC are special cases of the 
present 3-D formulation. The continuity of displacements and tractions, as well as the equilibrium 
equations, are imposed in the average (integral) sense as in the original HFGMC formulation. Each of the 
six sides (faces) of a subcell has an independent average displacement micro-variable vector which forms 
an energy-conjugate pair with the transformed average-traction vector. This allows generating symmetric 
stiffness matrices along with internal resisting vectors for the subcells which enhances the computational 
efficiency. The established new parametric 3-D HFGMC equations are formulated and solution 
implementations are addressed. Several applications for triply periodic 3-D composites are presented to 
demonstrate the general capability and varsity of the present parametric HFGMC method for refined 
micromechanical analysis by generating the spatial distributions of local stress fields. These applications 
include triply periodic composites with inclusions in the form of a cavity, spherical inclusion, ellipsoidal 
inclusion, discontinuous aligned short fiber. A 3-D repeating unit-cell for foam material composite is 
simulated. 

1.0 Introduction 

Early analytical three-dimensional micromechanical models have been concerned with 
homogenization methods able to generate the effective elastic properties of composites. A classical 
approach employs the Eshelby (1957) single inclusion elasticity solution in a homogenization scheme, e.g., 
the self-consistent scheme and its generalizations, Christensen (1979), and the Mori-Tanaka (Mori and 
Tanaka (1973)) methods. The latter approach has been widely used in the literature and extensions from 
single to double or nested inclusions have been reported, e.g., Hori and Nemat-Nasser (1993). Highlights of 
further extensions to compute the effective viscoelastic and inelastic composites has been given by Wang 
and Weng (1992) and Gavazzi and Lagoudas (1990), for example. Additional extensions to composites 
with inelastic matrix and progressive failure including debonding have been proposed by Ju and Tseng 
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(1996). The latter models have been generalized by Ju and Sun (2001) to develop a 3-D statistical 
micromechanical damage modeling framework for multiphase brittle composites with combinations of 
interacting microcracks and inclusions. Finally, combined damage 3-D micromechanics have been 
formulated and applied to functionally graded materials by Paulino et al. (2006). 

A different 3-D micromechanics approach is to use discrete based sub volumes where expansion of the 
displacement is carried out in order to capture the local fields. To this end, the classical displacement-based 
nonlinear FE method has been widely utilized, e.g., Levy and Papazian (1990). Jain and Ghosh (2008) and 
Jain and Ghosh (2009) used an extension from their 2-D to 3-D microstructural RVE model with Voronoi 
cell elements including damage. Multi-level variational formulation has been achieved to construct 3-D 
micromechanics with and without damage, e.g., Fish et al. (1999). Finally the previous micromodels and 
other simplified 3-D micromechanics have been used for local-global structural analysis. Towards that end, 
the lower-order GMC has been developed for 3-D composites. Thus, Haj-Ali et al. (2001) developed a 
nested 3-D micromodels for nonlinear multi-scale pultruded composite structures. Muliana and Kim 
(2007) employed a simplified 3-D micromodel using relatively small number of subcells for the 
viscoelastic multi-scale analysis of composite structures. The GMC has been utilized by Pahr and Arnold 
(2002) for the analysis of metal matrix composites with discontinuous fibers, by Bednarcyk and Arnold 
(2003) for the modeling of woven composites, and by Liu et al. (2011) for a multiscale analysis oftriaxially 
braided composites. The higher-order regular orthogonal 3-D-HFGMC has been employed by Bednarcyk 
et al. (2008) in their investigation of the foam insulation of the space shuttle. 

Refined micromechanics is a class of nonlinear micromechanical models capable not only to predict the 
instantaneous effective stiffness of a periodic heterogeneous medium, but also to accurately predict the 
local spatial distributions of the deformations within the repeating unit-cell (RUC). The latter requirement 
is needed in order to use such modeling for nonlinear material response along with potential localized 
damage, e.g., microcracking, fiber-matrix interface debonding, material defects and imperfections, among 
other modes. 

The E1FGMC method is a higher-order refined micromechanical theory and has evolved from its 
lower-order predecessors, the method of cells (MOC), Aboudi (1982) and Aboudi (1991), and the 
generalized method of cells (GMC), Paley and Aboudi (1992), respectively. The FIFGMC micromechanical 
model is well suited for detailed nonlinear analysis of multiphase composites, Aboudi (2004). Unlike 
general classical numerical methods such as finite element (FE) and finite-volume (FV), the E1FGMC is 
tailored and specialized to micromechanics of periodic composites. It is explicitly based on 
micromechanical variables needed to establish the elastic and inelastic concentration tensors of the phases 
in addition to the detailed local fields. The formulation of the HFGMC is performed using average 
equilibrium, traction and displacement continuity equations. This is because it aims to strike a balance 
between simplicity and accuracy without the need for excessive number of discrete sub-volumes (subcells). 
The computational affordability is especially important when HFGMC is integrated within a nonlinear 
analysis of composite structures including localized effects, often called multi-scale or local-global 
modeling. 

Several applications and implementations of the HFGMC have been performed to model different 
physical effects in periodic composites, such as coupled electro-magneto-thermo-elastic (Aboudi (2001)), 
viscoelastic-viscoplastic micromechanics Aboudi et al. (2002), Aboudi et al. (2003) and Aboudi (2005), 
bond damage of fiber, Bednarcyk et al. (2004) and Ryvkin and Aboudi (2007). A review of the method and 
application for HFGMC applied to smart materials can be found in Aboudi (2007) and Aboudi and Freed 
(2010). Damage evolution effects have been introduced within HFGMC by Haj-Ali and Aboudi (2009), 
Bednarcyk et al. (2010) and Aboudi (2011). These include cohesive growth between the subcells, cell 
extinction damage (CED), isotropic and anisotropic damage evolutions in unidirectional composites. 

The original HFGMC formulations have been performed in conjunction with orthogonal array of 
subcells used to depict the geometry of the phases (regular array). This limitation requires a relatively large 
number of subcells if the objective is to accurately capture refined geometrical features and the spatial 
variations of the local solution fields. Despite this limitation, the previous applications have demonstrated 
the effectiveness of the method to generate the local solution fields by using sufficiently large number of 
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rectangular subcells. It should be emphasized that few number of rectangular subcells is sufficient to 
generate the effective linear and nonlinear response for the overall composite with high accuracy. In fact, 
this is one of the advantages of using the HFGMC in a local-global analysis of composite structures with 
selective geometrical refinement for the microstructure, (Haj-Ali and Aboudi (2009)). 

A natural extension that overcomes the limitations of the regular orthogonal array is to employ a 
parametric mapping. To this end, Haj-Ali and Aboudi (2010) used linear geometrical mapping of the 
subcells to map the geometry of the phases of the composite. This linear and parametric geometric mapping 
can be applied for a unit-cell with general phase geometry using arbitrary quadrilateral cell shapes that are 
transformed to an auxiliary uniform square shape (natural coordinates). It is important to note that while the 
geometric mapping is linear, the subcell displacement expansion is of a full quadratic form (subparametric). 
In addition, the complete form of the Jacobian transformation has been employed without approximation 
which necessitates a numerical integration of the HFGMC equations. The utilization of a full quadratic 
expansion vector with its bilinear terms in a subcell required three additional equations which were chosen 
as the average moment of the equilibrium equations. It was emphasized by authors that this choice is not 
optimal and further research is needed to reach a proper form of the needed relations. 

In the present investigation, a generalization of the parametric HFGMC to 3-D multiphase composites 
with periodic microstructure is offered. To this end, general hexahedral subcells are introduced to model 
triply-periodic composites. Here too, the full quadratic displacement expansion is still used within the 
parametric HFGMC. However, we now show that the average displacement vectors at the edges (or faces) 
of the subcell can be expressed in terms of the coefficients of the expansion to form an independent subset 
of side-based variables. Furthermore, it is shown that corresponding transformed average -traction at the 
edges are conjugate pairs to the average displacements. In this way, two outcomes can be achieved. In the 
first, relations for the bilinear terms can be established with the quadratic coefficients of the displacements 
and thus no need for additional moment of equilibrium relations. In the second outcome, the energy 
conjugate pair result in symmetric linearized form of the overall HFGMC equations with obvious 
computational advantages. The present new 3-D parametric HFGMC generalization allows the analysis of 
composites with general shape of inclusions. The 2-D parametric HFGMC can be obtained as a special case 
of the present formulation. 

This paper is organized as follows. In Section 2.0, the general theory to 3-D parametric HFGMC is 
presented including the 2-D-HFGMC as a special case. In Section 3.0, computational implementation 
aspects are described and discussed. Wide range of applications are given in Section 4.0, followed by 
concluding remarks in Section 5.0. 

2.0 General Formulation 

The HFGMC micromechanical method is presently extended for general parametric formulation 
suitable for three-dimensional (3-D) analysis of multiphase composite materials. Figure 1 schematically 
illustrates a triply-periodic multiphase material system having a global coordinates (jci, x 2 , x 3 ). The 
repeating unit cell (RUC) of this periodic composite can be identified and described by using the coordinate 
system (y u y 2 ,y^)- The goals of the HFGMC micromechanical method, based on the homogenization 
technique for periodic composites, are to predict the overall effective properties as well as the spatial local 
deformation fields. In the HFGMC framework, the RUC is divided into array of cells, often denoted as 
subcells. In the present section, this method is extended to analyze an RUC of a triply-periodic composite 
with irregular array of generally shaped hexahedral subcells. Figure 1 shows a schematic RUC domain 
which is discretized into a general assembly of hexahedral subcells to represent the different phase 
geometries. Figure 2 shows a general hexahedral subcell isolated and shown in its physical coordinates 
Oi, y 2 , v 3 ). This subcell is mapped to a uniform parametric coordinate system (r, s, t) using the classical 
linear transformation 
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z =1,2,3 


( 1 ) 


8 

y i (r,s,t) = ^ H k (r,s,t)y k i, 
k = 1 

where y h i = 1, 2, 3, are the coordinates of a general point within the subcell ((3) mapped from the parent 
parametric coordinates to the physical RUC coordinates, and 

H x = 1(1 - r)( 1 - s) ( 1 - t), H 5 = 1(1 - r)( 1 - ,s’)(l + f) 

O O 

H 2 = 1(1 + r)(l - ,S')(1 - 0, H 6 = 1(1 + r)(l - .S’)(l + 0 

O O . 

H 3 = i(l + r)(l + s)( 1 - 1), H 7 = |(1 + r)(l + s)( 1 + 1) 

O O 

// 4 = |(1 - r)(l + ,s')(l - 1), H 8 = |(1 - r)(l + s)( 1 + 0 

O O 

and Vki are the comer coordinates of the subcell. 



Figure 1 . — Schematic illustration of a triply periodic array in 
the global (x-i, X2, X3) space of multiphase composite media 
with its repeating unit-cell (RUC) having hexahedral 
subcells, defined with respect to its (y-i, /2, 73) local 
coordinate system. 
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Figure 2. — A hexahedral shaped subcell in its physical and 
natural coordinates. Numbering for the faces and corners are 
shown. 


In the present 3-D HFGMC formulation, the complete quadratic form of the displacement expansion in 
the subcell, is given by 

u = £ • x + W 00 o + W 100 r + WqjqS + W 001 f + Wj l0 rs + W 101 rt + W 01 ^st 
+ j W 200 (3r 2 - 1)+ j W 020 (3 s 2 - 1)+ jW 002 (3t 2 - 1) ^ 

where the applied remote global-scale displacement field is given by u° = e () • x. This expansion is a direct 
extension of the full quadratic expansion which has been used by Haj-Ali and Aboudi (2010) in the case of 
the parametric HFGMC for doubly-periodic composites. However, in the present parametric formulation, 
the mixed bilinear terms of this expansion are shown to be dependent variables. This resolves the open 
question that has been raised by Haj-Ali and Aboudi (2010) on the optimal method that should be used to 
determine these terms. Previously, the bilinear terms have treated as independent variables and 
consequently additional moment of equilibrium equations have been introduce to determine them. To this 
end, let us define the average displacement vector on the six faces of the hexahedral subcell: 

u ( M=_Lf u (| 3 ) (y)cM yt =-jf f u^\^ k ,r\ k )dt lk dr\ k , k = 1 , 2,.., 6 ( 4 ) 

A 4 J- 1 J -1 
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where (£#, r)*) are the surface parametric integration variables of the k - th edge or face of the subcell denoted 
by (Pi). The above integration provides the six vectors of the surface-average displacements expressed in 
terms of the expansion coefficients W ( /„ m) as follows: 




u°,(M 


w,~ 


' 0 

-1 

0 

0 

1 

0 

f 


W 100 ' 



„0,(P 2 ) 


W 2 


1 

0 

0 

1 

0 

0 

1 


W 010 

u<fc> 


u«,(p 3 ) 


W 3 


0 

1 

0 

0 

1 

0 

1 


Wool 


» — < 

„0,(P 4 ) 

> = « 

W 4 

■ = 

-1 

0 

0 

1 

0 

0 

1 


W 20 o 

U<0 5 > 


u°.<p 5 ) 


W 5 


0 

0 

-1 

0 

0 

1 

1 


W 0 20 



u°’(P‘) 


W 6 


0 

0 

1 

0 

0 

1 

1 


W 0 02 

w 000 . 


0 


Wo. 


0 

0 

0 

0 

0 

0 

1 


W 000 . 


The inverse relation of Equation (5) is given by: 


w 100 ' 


' 0 

1 

0 

-1 

0 

0 

0 " 


wr 

W 010 


-1 

0 

1 

0 

0 

0 

0 


w 2 

W 001 

_ 1 
’~2 

0 

0 

0 

0 

-1 

1 

0 


w 3 

W 2 oo 

0 

1 

0 

1 

0 

0 

-2 


w 4 

W 0 20 


1 

0 

1 

0 

0 

0 

-2 


W 5 

W 0 02 


0 

0 

0 

0 

1 

1 

-2 


W 6 

W 000 . 


0 

0 

0 

0 

0 

0 

2 


w o. 


( 5 ) 


(6) 


It can be readily observed that the bilinear coefficients Who, Wioi, Won, do not contribute to the average 
displacements. Thus it is possible to assume that these are dependent coefficients if our aim is to use the 
average displacements as the primary independent variables. Hence, two possible solutions can be 
identified. The first trivial case is to assume that the bilinear terms to be zero: 


W 110 =w 011 =w 101 =0 ( 7 ) 

This solution generates a truncated quadratic expansion that may affect the spatial distribution of the elastic 
fields. However, the average strain in the subcell, directly responsible for the overall elastic effective 
stiffness of the RUC, is not affected due to the vanishing integration of the bilinear terms in the hexahedral 
domain. 

Alternatively, a non-trivial solution can be obtained while retaining the complete quadratic expansion 
(3), maintaining terms symmetry and frame indifference of this polynomial form. Here 

W, l0 = ^-(W 200 + W 020 ) = 1(W 3 + W 4 + W 5 + W 6 - 4 W 0 ) 

Won = |(W 020 + W 002 ) = ±(Wj + W 2 + W 3 + W 5 - 4W 0 ) (8) 

W 10 i = ^(W 200 + W 002 ) = i( Wl + W 2 + W 4 + W 6 - 4 W 0 ) 

It can be verified that the average displacement relations given by Equation (5) is maintained with the trivial 
and the non-trivial forms given for the bilinear terms. Equations (7) and (8), respectively. It should be noted 
that the trivial case is a direct result of the original HFGMC formulation, e.g., Aboudi (2004), for 
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orthogonal and regular array of subcells where both physical and natural coordinates coincide and the 
Jacobian of the mapping is constant. 

As a result of the above derivations it is possible to rewrite 
following compact form: 

u = e°-x + W 0 + ^(W 4 - W 6 )r + i(W 5 - 

+ ^-(W 4 + W 6 - 2W 0 )(3r 2 + rs + rt - 1) 

+ i(W 5 + W 3 - 2W 0 )(3s 2 + rs + st - 1) 

+ ^-(Wi + W 2 - 2W 0 )(3 1 2 + rt + st - 1) 

It is worth mentioning that in the special case of 2-D parametric HFGMC for doubly-periodic composites, 
the above expansion takes the form: 

u = 8° • x + W 0 + ^-(w 2 - W 4 )r + i(w 3 - wjj 

+ i(W 2 +W 4 -2W 0 )(3r 2 +rs-l) (10) 

+ ^-(W, + W 3 - 2W 0 )(3s 2 + rs - 1) 

where W* follows the labeling of the quadrilateral subcell of Haj-Ali and Aboudi (2010). In this case, the 
bilinear term coefficients Wn takes the following trivial and non-trivial forms, respectively: 

W n =0, Wj j = ~(w 20 + W 02 ) (11) 

As mentioned above, this answers the open question raised by the authors on the optimal way of handling 
the bilinear terms in the 2-D case. 

In order to establish the expressions for the strain components, the displacement gradients with respect 
to the physical coordinates y are needed. The standard Jacobian of the linear transformation (1) is given by: 


the displacement expansion (3) in the 


W 3 > + |(W 2 -Wj> 


(9) 


j _ dOTJTDb) 

d(r,s,t) 


dy, 

dy 2 

8y 2 

dr 

dr 

dr 

dy, 

8y 2 

8v } 

8s 

8s 

8s 

dy, 

dy 2 

dy 2 

dt 

dt 

dt 


( 12 ) 


Its inverse is denoted by J =J 1 = d(r,s,t)l d(y x ,y 2 ,y^) and used to relate the displacement derivatives 
as 


r 

(P) 



r i 

u i 




U ,r 



(P) 


U 2 

' — 

J 

< 

U .s 

u ,3 




U ' 


The right-hand -side can be easily established by using Equation (3) yielding 


(13) 
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( 14 ) 


= W 100 + W 110 s + W 101 t + 3W 200 r 
< = W 010 + Wj 10 r + W 011 t + 3 W 020 s 
= W 001 + W 101 r + W 011 s + 3W 002 t 

However, in order to express these gradients in terms of the average displacement micro-variables of the six 
hexahedral faces, we first define the following subcell vector with its components lumped in the following 
order 


w (w ’ r = {W u ,W 2i ,W 3i ,W 4i ,W 5i ,W 6i ,W 0i } m 1 = 1,2, 3 

Next, Equations (6) and (8) are used for the non-trivial case to obtain, after some manipulations, the 
displacement gradients with respect to the parametric variables as follows 


(P) 


= [M] (P) 


W, 

W 2 

w 3 

w 4 

W 5 

W 6 

Wn 


(P) 


where the matrix [M] I|S| is given by 


s -6 r - 4(5 + 1 ) 
r + t - 6s - 4(r + 1) 
— 1 + r + 3t 1 + s + 3t —6 1 - 4(r + s ) 

Similarly, Equations (6) and (7) provide the gradients for the trivial case as 


t 1 + 3 r + t s + t — 1 + 3r + s s + t 

- 1 + 3s + 1 t 1 + r + 3s r + t r 

r + s r + s r s 


U ,r 

(P) 

1 

o 

0 

1 + 3 r 

0 

- 1 + 3 r 

0 

0 

- 6 r 

U ,v 

-l + 3s 

0 

1 + 3s 

0 

0 

0 

- 6s) 

U ,1 

Z 

0 

0 

0 

0 

-1 + 3* 

1 + 3* 

- 6 1 


w, 

W 2 

W 3 

w 4 

W 5 

W 6 

w n 


(P) 


The spatial form of the strains in a given cell ((3 ) are defined by 


,(P), 


’(y(^,s)) = ej +j- 
It is possible to bring the strains into a general matrix form 


1 


dU: Su . 


yfyj tyi 


(15) 


(16) 


(17) 


(18) 


(19) 
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where the strain vector notation s is defined by 

£ — { £ 11j £ 22» 8 33’^ 8 23’^ 8 13’^ £ 12} (^1) 

and A (P) is the matrix that relates the strain to the displacement microvariables. The size of the vector W (,1) is 
21 which includes all the microvariables of subcell ((3) as shown in Equation (15). Using Equations (17) and 
(18), respectively, along with Equations (13) and (18), it is possible to derive the strain-displacement 
matrix, A (P) . The non-zero components of this matrix are listed in the Appendix for the trivial and 
non-trivial bilinear terms. 

Next, we define average traction vector T (Pa) on the six sides of the hexahedral subcell defined by 


Tdh) — _j_ | 


[ CT (P) (y)-n (Pi) i/5 A , f [ a m (r,s,t)-n ih) d^ k dr\ k , k = 1,2,.., 6 (22) 

JSb 4 J— U— l 


where n (Pi ^ is the unit nomial vector to the k - th side of subcell (|3). The stress field, a <P) , is given by 

CT ( P) _ c ( P* ; g ( P ) (23) 

assuming a linear material behavior of the subcell, with C <P) being the stiffness of the material in the subcell. 
Since the geometry is interpolated using linear mapping, the normal vector n (p * ) to each of the six faces is 
constant. Employing the established spatial strains in the cell to obtain the corresponding stresses. The latter 
are used in the expression for the average tractions, Equation (22). This provides 


if f 

4 J-iJ-i 


r p(Pi) _ J_ I " 1 I " 1 jy(Pt)^(P) 


with, 


jy(Pt) 


and 


] ]d^k dr \ 


8 ° + A (P) W (P - An ~ lvlPilrtP) 


N' PA 'C' 


e° + A (Pi) W (p) 


(24) 


n x 0 0 0 « 3 n 2 

0 « 2 0 « 3 0 

0 0 « 3 n 2 «i 0 


(PU 


(25) 




(26) 


As mentioned above, due to the applied linear mapping, the matrix is constant. 

Following the original HFGMC formulation, we impose the displacements and tractions continuity 
between the subcells. These conditions are enforced in an average integral sense. In addition, periodicity 
conditions are imposed between the boundary subcells of the RUC by requiring that the displacements and 
tractions be continuous at these mirrored interfaces. The third major requirement of the FIFGMC method is 
the intra subcell equilibrium applied in a volumetric average form. In the following, the above transformed 
expressions for the average displacements and tractions, (4) and (24), respectively can readily be used to 
impose the FIFGMC equations. The displacements and tractions continuity can be written as 
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( 27 ) 


y(Pt) _ ^jyy(P) _ y(Y m ) _ \y (y) = r p('l'm) 

where (P^) denotes the k - th interface (side) of subcell ((3) and (y m ) is neighboring /77-th interface side of 
subcell (y). The displacement and traction periodicity conditions are imposed as in Equation (25), but with 
cell (P) and (y) located on opposite sides of the RUC. 

The equilibrium equations for each subcell are imposed in an average sense in conjunction with 
Divergence theorem in order to utilize the derived expressions for the average tractions as follows. 

6 6 

J V-a</F = J’ a-ndS = N (Pt) a (P) dS t =j^*T (Pi) =o (28) 

^ / — 1 7 1 . 


where V and S are the volume and surface of the (P)-subcell, respectively, and A k is the area of the k - th side. 

The total number of unknown microvariables in the RUC is 21 N c where N c is the total number of 
subcells every one of which has 21 microvariables, see Equation (5). The number of continuity and 
periodicity equations for displacements is 3*3N C (only three interfaces for each subcell provide 
independent relations), see the first equality in Equation (27). Similarly, the number of traction continuity 
and periodicity equations is 3x3 N c as shown by the second equality in Equation (27). The average 
equilibrium relations. Equation (28), provide additional 3N C equations. 

The above formulation completes the full theoretical derivation of the general three-dimensional 
parametric HFGMC micromechanical method. However, it is interesting to introduce an internal force 
vector I (P) that form energy conjugate pairs with the average displacements of each side of the subcell, W (P) . 
The internal force vectors are directly related to the average traction defined above. These force resisting 
vectors are defined using external and internal work balance: 

W (P)^I<P) s ( e (P)T _ g °T ) o m dv = ^ w (P) - r A (P) ' r a (P: W ( 29 ) 


and since W is arbitrary, the expression for the internal resisting vector is 

I(P) = f A (P)J o (P Vf 

Jv 


(30) 


The above work balance is known in the mechanics of heterogeneous media as the principal of average 
virtual work, e.g., Christensen (1979). For the linear material case, 


r(P) _ 


4 


,(P)T, 


Cdv 


| £ °,(P) + 


[[ A^CA^rfv 


W 


(P) 


(31) 


We further recognize that the first six vectors within I (|3) are the average forces, , on each face or side 
of the hexahedral subcell. The use of the above form of the average internal resisting vector, I ( ^ , instead 
of the average traction, , in the traction continuity relation and average equilibrium, Equations (27), 
(28), respectively, provide enhanced computational advantage mainly expressed by using the symmetric 
matrices in Equation (31). 


3.0 Implementation and Computational Aspects 

The assembly of the derived 21 *N C HFGMC governing equations of the RUC are solved to obtain the 
unknown 2 1 microvariables for each cell. In general, equations can symbolically be grouped into three parts 
in the form 
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( 32 ) 


1 

i 

i 


0 


0 

4 +) +i^ 

{w} + 

D^ +) 

II 

O 

0 

i 

M 
1 


1 

M 

d 

1 


0 


where the first part (row) represents the average continuity and periodicity of the displacements, the second 
part represents the continuity of the tractions in a similar fashion. The third part represents the equilibrium 
equations for all the cells. The above system of equations is solved for a given externally applied strain 8° to 
obtain the cell microvariables. Furthermore, the local spatial strains can be readily obtained using 
Equation (20), which leads to the spatial distribution of the stresses in the RUC. These can be used, for 
example, in a progressive analysis to determine possible damage in the constituents as well as interfacial 
debonding, see Haj-Ali and Aboudi (2010) for the parametric 2-D-HFGMC case. 

Since the first part of these equations are homogeneous and each row is composed of only two terms, 
expressing interface average displacement continuity becomes a straight forward manner. A pre -analysis 
overall RUC equivalence of matching subcells' faces can be performed similar to nodal equivalence earned 
out in the FE method. This can be done at the global level of the code or even using simple preprocessing 
subroutines. Furthermore, the equilibrium equation for each subcell. Equation (28), expressed in either 
T (Pt ) or j(PD ^ can [-, c usec j j n conjunction with a static condensation of the internal dependent 
displacement microvariables, W 0 <P) . This leaves the six-side average displacements as the only set of 
independent variables for each subcell, which results in computational saving. 

The effective elastic properties of the multiphase composite can be obtained from the cell strain 
concentration tensors G ( |l! which relates the average strain in the subcell to the externally applied strain. The 
latter is established by considering the average strain of the cell having the form 

= bj ' r s(y)JV = e» + Tj (v,„ + n V r )dV (33) 

By using the Divergence theorem and the linear mapping. Equation (30) takes the form 

1 i 6 

8 (P) = e° + — | (u®n + n(8)uV5 = s 0 + — ^4(u®n + n®u) (pA) =G (P) : 8° (34) 

2V s 2V k=l 

It is observed from Equation (34) that the concentration tensor G l|; " depends on the average displacement 
vectors on the surfaces of the cell. Those can be computed only after solving the entire FIFGMC system of 
equations for the RUC as discussed above. Once the concentration tensors G (P) have been obtained, the 
effective elastic stiffness tensor C is evaluated in the form 

C*=f>< P) C (P) G ,P) (35) 

p=i 

where v p = 47 V RUC , and V RUC being the total volume of the RUC. 

It is also possible to use the intact overall system of equations. Equation (34), without side-based 
condensation between the displacement microvariable. This approach allows the formulation of a new 
micromechanical damage approach, termed cell-extinction-damage (CED). As has been demonstrated by 
the 2-D-HFGMC, Haj-Ali (2009) and Haj-Ali and Aboudi (2010), this cohesive micromechanical 
approach has the ability to perform subcell degradation along with general traction-separation between the 
subcells. In addition, the structure of exclusive interface-based equations unlike nodal-based FE, makes it 
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easier to remove the subcells from the problem and/or generate average -based traction-free surfaces for 
crack propagation within the RUC without introducing computational instabilities. 

As mentioned, linear geometric mapping and quadratic displacement interpolation (subparametric 
formulation) has been carried out in the previous section. This formulation has been done in order to 
simplify the numerical integration. Having said that, extension to full quadratic geometric mapping 
(isoparametric formulation) is straightforward. In both of these cases, numerical integration of the 
equations is needed since the Jacobian matrix is not constant. However, a computational implementation of 
the 2-D parametric HFGMC has been reported in the literature under the name “FVDAM”, e.g., Bansal and 
Pindera (2006). This latter 2-D method assumes a volume-based constant Jacobian of the linear mapping 
for the geometry of the subcells. Their assumption imposes severe limitation and restriction on the 
parametric HFGMC method as the shape of the subcells need to be close to rectangular or parallelograms in 
order to maintain accurate depiction of the geometry. The use of an a priori constant Jacobian amounts to a 
one center-point quadrature and may lead to large solution errors especially in the displacement gradients. 
Furthermore, the assumption of a constant Jacobean is not applicable in the general 3-D HFGMC that has 
full quadratic displacements. Thus, the proposed 3-D parametric HFGMC has been formulated in a direct 
and general manner, i.e., retaining the full quadratic expansion of the displacement together with the 
complete Jacobian. The level of computational accuracy is determined by the number of integration points 
specified during the numerical integration of the governing equations. 

It is important to draw the distinctions between the proposed parametric HFGMC and the classical 
displacement-based finite-element (FE) method, where a common misconception has been to link the 
HFGMC to FE. In the FE formulation the displacement continuity between two adjacent and connected 
elements is satisfied in a pointwise manner by sharing the same nodal degrees of freedom at the sides of an 
element. However, the HFGMC quadratic displacement expansion is nonconforming and displacement 
continuity is satisfied in an average sense between adjacent cells. The latter is an approximation that allows 
the HFGMC to explicitly use additional stress variables in the formulation and directly apply average 
traction continuity. The advantage of this approach is that both traction and displacement interface 
continuity between different materials can be directly employed. Further, cohesive damage modeling can 
be easily employed as mentioned before for the 2-D case. Unlike FE, the HFGMC is a specialized 
micromechanical modeling framework in which the remote fields are directly tied to the micromechanical 
formulation (e.g., the remote average strain is directly used in the HFGMC formulation). On the other hand, 
the periodicity in the FE formulation is imposed only through displacement boundary conditions, while in 
the HFGMC, the periodic boundary conditions are imposed directly using both the average displacement 
and tractions. The latter differences are illustrated in the next application section where examples for RUCs 
with both traction and displacement periodic conditions are compared to FE RUCs modeled with applied 
far-field nodal displacements. Having said that, it is possible to draw an analogy between the FE and 
HFGMC method by considering the average displacement micro variables as face or side-based generalized 
displacements variables. This analogy allows the use of the average virtual work principal to define the 
internal force vector and its associate symmetric matrices which are useful for the solution of the overall 
system of equations. 

4.0 Applications 

Several applications are performed to show the effectiveness of the parametric HFGMC method in 
capturing the spatial stress field distributions within the phases of triply periodic composites. To this end, 
RUCs are generated with coarse and refined resolutions using hexahedral subcells. Selected stress fields are 
compared with solutions of RUCs using displacement based FE models. In order to compare the parametric 
HFGMC and the FE solutions, the full periodic conditions between the boundary surfaces of the RUC, i.e., 
both traction and displacement continuity on mirrored boundary surfaces, were not applied. It is important 
to mention that the full periodic conditions are an integral part of the HFGMC formulation, however, full 
periodicity is not a straightforward using the displacement-based FE. Thus, the boundary conditions 
imposed in all compared cases, unless stated otherwise, are: 


NASA/CR— 2012-217715 


12 



U°(5) = 8° -x 


(36) 


where (S) here is the outside periodic surface of the RUC, e° is the applied average strain and x are the 
coordinates of the surface points. In the case of the FE-RUC, constrains were written for the boundary 
nodes that link their degrees of freedom to the applied average strain. However, for the case of the RUC 
modeled with the HFGMC, the average displacement microvariables on the subcells with boundary faces 
were set to zero. Select cases for HFGMC-RUCs with full periodic boundary conditions are presented in 
order to show the different spatial distribution of the local fields that the can be achieved using the HFGMC. 

The first RUC model is for a medium with repeated spherical cavities, thus forming a porous material. 
Figure 3 shows the geometry of this problem in the form of a half RUC cut at the mid-plane and another 
with a quarter portion removed. It should be mentioned that that while it is possible to apply symmetric 
conditions, the full geometry is simulated. The dimensions of the RUC were taken as D = H = L=\ and the 
radius of the spherical cavity was R = 0.4. The material Young's modulus and Poisson's ratio elastic 
properties are E = 4.8 GPa, v = 0.34. The composite is subjected to a remote uniaxial strain loading of 

= 0.01 . The number of hexahedral subcells used to map the geometry was N c = 3584. This discretization 
of the RUC is not considered as a refined one compared to the number of subcells needed to capture the 
local spherical shape curvature. This was the choice in order to demonstrate the ability of the HFGMC to 
accurately solve for the local fields despite the use of average traction and displacements. Figure 4(a) and 
(b) show the an stress distributions for the HFGMC and the FE models, respectively. In these two cases, 
only the displacement-based periodic conditions were enforced. It is clear that the stress distribution from 
the HFGMC solution is quite identical to the one obtained from the FE model despite the use of relatively 
smaller number of independent solution variables. As mentioned, each subcell has 2 1 N c total independent 
average displacement micro variables in the case of the HFGMC among of which many dependent 
displacements are condensed due to subcells inter-continuity at the global system of equations as discussed 
in the implementations section. Therefore, the element that was selected for the FE RUC cases is the 
8 -nodes linear brick which has 8><3 *N C = 24 N c number of variables having both RUCs close in their 
problem size. Figure 4(c) shows the stress distribution generated by the HFGMC, however, the frill periodic 
conditions were applied in this case. It is evident that the local spatial distribution is affected at the 
boundaries of the RUCs if accurate periodicity needs to be captured. In the case where the effective elastic 
properties of the entire medium are sought, applying frill or partial periodic conditions may not yield 
pronounced differences. Our aim in the parametric HFGMC is to provide an analysis tool capable of 
accurately predicting the local fields. 

The next application deals with a composite with spherical inclusions, thus forming a particulate 
composite. Figure 5 shows the discretized RUC with D = H = L = 1 and the radius of the spherical inclusion 
is R = 0.4, which corresponds to inclusion volume fraction of 26.8 percent. An effort has been made to use 
a refined mesh at the interface of the spherical inclusion. The elastic properties of the inclusion and the 
matrix are: E = 380 GPa, v = 0. 1 and E = 4.8 GPa, v = 0.34, respectively. The far field applied strain is 
833 = 0.01 . The number of hexahedral subcells used for both HFGMC and FE is N c = 74088. Figure 6 (a) 
and (b) show the Mises stress distributions as predicted by the HFGMC and the FE, respectively, in the 
direction of the applied strain. Here, excellent agreement between the two solutions is exhibited. As a 
continuation of this problem, we consider an ellipsoidal inclusion, the mesh of which is shown in Figure 7. 
It illustrates a half RUC with removed inclusion and the ellipsoidal inclusion. A coarser mesh is used for the 
latter in order to minimize the size of the problem: N c = 54872. The dimensions of the RUC are D = H = L 
=1 and the major and two minor axes of the ellipsoid are 0.8, 0.4 and 0.4 (spheroidal shape), respectively. 
The same elastic properties for the inclusion and matrix as well as the magnitude of the applied far field are 
chosen as in the previous problem. Figure 8 displays the resulting stress a 33 distribution in the major-axis 
direction of the ellipsoid for both the HFGMC and FE solutions and for the full RUC with quarter portion 
removed. Again, a very good agreement exists. 
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Figure 3. — Parametric HFGMC meshes for an RUC with spherical cavity. 



Figure 4. — Stress distributions in the RUC with a spherical cavity, (a) HFGMC results (periodic displacements), 
(b) FE results (periodic displacements), and (c) HFGMC results (periodic displacements and tractions). 
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Figure 5. — Parametric HFGMC 3-D meshes for an RUC with spherical 
inclusion using hexahedral subcells. 



Figure 6. — Stress distributions in the RUC with a spherical inclusion, (a) HFGMC 
results, (b) FE results. 
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Figure 7. — Parametric HFGMC 3-D meshes for an RUC with 
spheroidal inclusion using hexahedral subcells. 



Figure 8. — Stress distributions in the RUC with a spheroidal inclusion, (a) HFGMC results, 

(b) FE results. 

The next studied problem analyzes a composite with a system of reinforcements in the form of 
discontinuous aligned circular fibers (short-fiber composite). Figure 9 displays the refined mesh of the 
RUC together with the short fiber itself and its cross-section. The ratios of the dimensions of the RUC are 
0 = 5 and H = L= 1 . The axial length of the fiber is 4.53 and its diameter is equal to 0.4, resulting in a fiber 
volume fraction of 1 1.4 percent. The RUC is subjected to remote axial strain = 0.01 . The values of the 
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fiber and matrix elastic properties are the same as those of the inclusion and matrix properties in the 
previous problem. The size of the problem in this case is N c =17984. The resulting Mises stress distribution 
within the RUC as predicted by the HFGMC and the FE are shown in Figure 10(a) and (b), respectively. It 
can be observed that the FIFGMC is capable in capturing the local mechanism of stress-transfer at the 
discontinuous edge of the fiber to the matrix. Both solutions are very close as shown. 

The final example considers a periodic foam-like material with a relatively complex geometry of the 
microstructure. Figure 1 1 shows the selected geometry of the RUC. It is evident that this is a relatively 
complex geometry due to the spatial curvature and shape of the cuts. Accurate modeling of the material in this 
case requires a refined 3-D mesh able to represent the local geometry. The choice of the RUC geometry has 
been inspired from scaffolds used in tissue engineering and orthopedic micro-implants. Its dimensions are D = 
H = L= \ and was constructed by removing material enclosed by three orthogonal cylinders, the diameter of 
each is 0.8. It is worth observing the local geometry where the removed cylinders intersect. The material 
elastic properties are E = 4.8 GPa, v = 0.34. The total number of subcells is N c =34560 and the remote axial 
strain is = 0.01 . Figure 12(a) and (b) exhibit the Mises stress distribution in the 3-D FIFGMC compared 
with the FE solution, respectively. Excellent agreement exists especially at the inner comers. It is important to 
mention for this case, that an explicit average traction-free conditions are imposed by the FIFGMC on the 
inner surfaces. Finally, the Mises stress variations are shown in Figure 13 for the same RUC but when full 
periodicity conditions are applied between the RUCs by the FIFGMC model. It is clearly observed that the 
fomi of the stress distribution is impacted by imposing the full micromechanical boundary conditions. This 
demonstrates the micromechanical capability embedded in the formulation of the FIFGMC. By introducing 
inelastic effects, it would be possible to extend the 3-D parametric FIFGMC to analyze metallic foams. 



Figure 9. — Parametric HFGMC 3-D meshes for an RUC, with discontinuous 


aligned fiber having a circular cross-section, using hexahedral subcells. 
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Figure 10. — Stress distributions in the RUC with the discontinuous fiber, (a) HFGMC 
results, (b) FE results. 
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Figure 1 1 . — Parametric HFGMC mesh for an RUC of a foam material. 



Figure 12. — Stress distributions in the RUC of the foam material, (a) HFGMC results (periodic displacements), 
(b) FE results (periodic displacements). 


NASA/CR— 2012-217715 


19 



Figure 13. — Stress distributions in the RUC of the foam 
material showing HFGMC results with full periodic 
conditions (displacements and tractions). 


5.0 Conclusions 

A new 3-D parametric micromechanical formulation of the HFGMC method is presented using linear 
geometric mapping with hexahedral subcells. A complete quadratic expansion of the displacement is used 
for the subcells where its u nk nown coefficients are in the form of average displacements of the six sides of 
the subcells in addition to an internal coefficient that can be statically condensed using the integral-based 
equilibrium equations. The independent average displacement micro-variable vector for each side of the 
subcell forms an energy-conjugate pair with the transformed average-traction vector. Trivial and non-trivial 
relations for the coefficients of the bilinear terms are found and their corresponding strain-displacement 
matrices are derived without the need for additional equations in the form of moments of the equilibrium. 
Thus, the open question previously posed by the authors on the optimal equations for the bilinear 
coefficients is resolved in the current study. The previous 2-D parametric and orthogonal HFGMC 
formulations are special cases of the present 3-D formulation. Solution approaches to the parametric 
HFGMC equations are also discussed and numerical integration is needed due to nonconstant Jacobian of 
the linear mapping. Wide range of applications are presented to verify the 3-D parametric HFGMC as 
compared with finite-element RUC models. The proposed 3-D parametric HFGMC is general and well 
suited for nonlinear material and finite-deformations previously developed in the framework of the 
orthogonal HFGMC. Previously developed cohesive damage modeling can be integrated within the new 
formulation. These modeling features can be integrated within the proposed 3-D parametric 
micromechanics along with local-global multi-scale analysis of composite structures. 
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Appendix 

In this Appendix, the non-zero components of the strain micro-displacements matrices are listed first 
for the quadratic polynomial with zero bilinear terms (trivial case) and second for the constrained 
symmetric bilinear term coefficients (non-trivial case). The terms for first matrix are 

A u = -1/2 Jj 3 + 3/27 1>3 / 

4 — 1 / 27 1 3 +3/2A 3 / 

A l7 = -1/2 J X 2 +3/2 J l2 s 

T uo =1/2 Jj j +3/2 J l x r 

A ll3 = 1/2 J l2 + 3/2 J l2 s 

A lA6 =-1/27]] +3/2 J lA r 

A\ ,19 = _ 3 J\ \Y - 3 J\ 2 S ~ 3 A, 3 

t A, 2 =-l/27 2i3 +3/27 2 ,3t 

A, 5 = 1 / 2 J 2 ,3 + 31 2J 2<3 t 

A ,8 = -l/27 22 + 3/27 2j2 s 

A 2AX = 1/27 2i1 +3/2 J 2l r 

T 2 14 = 1 / 2 7 2 7 + 3/2 J 22 s 

A , 17 = -1/2 7 21 +3/2 J 2l r 

A~t 2 q — 3 7 2 ]7* 37") 2 5 3J 2 3 1 

A 3 3 = — 1 / 2 J i + 3/ 2 J 3 3 t 

A 3 g = 1/27, 3 + 3 / 2 7 3 3 / 

A 39 = -1/2 J 32 + 3/2 J 3t 2 s 

A , 12 = 1/2731 + 3/2J 3A r 

A 3l5 = 1/2 J 32 + 3/2 J 3 2 s 

A, is = — 1 / 2 A 3 ] +3/2J 3A r 

A 3 21 = -3 J 3A r - 3 J 3 2 s -3 J 3 3 t A 42 = -l/2J 33 +3/2 J 3 3 t 

A 4 ' 3 = -1/2 J 23 +3/2 J 2 3 t 

A 4 5 = 1/2 J 3 3 + 3/2 J 3 ' 3 t 

A 4 g = 1/2 J 1 3 +3/2 J 2 3 t 

A 48 = -l/2J 32 + 3 / 2 7 3>2 s 

A 4g = -1/2 J 22 + 3/2 J 2 _ 2 s 

7 4 ,h =1/27 3 j +3/2J 3A r 

A, 12 =1/2 J 21+ 3/ 2 J 2 \Y 

A 4A4 =1/27 3 i 2 +3/27 3j2 s 

A 415 = 1/2 J 2>2 +3/2 J 2 ' 2 s 

A , 17 =-1/27 3 ,+ 3/27 , p* 

A 4Ag = -1/2 J 2A +3/2 J 2A r 

A , 20 = -3 7 3 p' — 3 7 3j2 5 — 37 3 3 / 

2 i — ^ J 2\^ 3 ^2 2 *^ ^*^2 3^ -^5 i — 1 / 2 «y 3 3 3 / ^ A ^ 

A 5 ' 3 = -\/2J l3 +3/2J l3 t 

T, 4 — 1/2733 + 3 / 2733 / 

A 3 6 = 1/27, 1 + 3 / 2 7, ,/ 

T, ~j = — 1/2 7 3 7 + 3/2 J 3 2 s 

A 59 = -1/27, 2 +3/2 J l2 s 

A s ]q =1/27 3j1 + 3/27, p* 

A 512 = 1/2 7, , +3/2 J lA r 

A, 13 =1/27 3j2 +3/27,25 

A 5A5 = 1/2 J l2 + 3/2 J l2 s 

T 5 16 =-1/27 3j1 +3/27, p* 

As, is = -1/2 j + 3/2 J 1A r 

T, ]g = 373 , 7 * 37 3 2 5 3 A 3 3 t 

A, 21 ~ “3 A,i r — 3J X 2 s — 3 J l 3 

t A,i=-l/2i 2 ,3 +3/27 2>3 / 

Afi 2 = -1/2 J, 1 +3/2 7, 3 t 

^6 4 =l/2 7 2 ,+3/2 7 2 3 / 

Afr 3 = 1/2 A 1 3 +3/2 J \ 3 t 

A , 7 =-l/27 2 2 +3/27 2 2 5 

A 6% = -1/2 J X1 +3/2J hl s 

A,io = 1/ 27 2j ] + 3/27 21 7* 

A 6U = 1/2 7, , +3/2 J lA r 

T 6 13 =1/27, 2 + 3/2 J 22 s 

a 6A 4 = 1/2 y 1>2 + 3/2 J 12 5 

A,i6 = -1/ 27 2 ] + 3 / 2 J 2A r 

T 6 17 = -1/2 A j + 3/2 J lA r 

A,i 9 = T 2j ] 7* - 3 7 2 2 s - 3 J 2 ^t 

A , 20 = -3Tj p' - 3 7, 2 s - 3 

t 
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The non-zero terms for the second matrix are 


A u = l/4J l3 r +l/4J 13 s+(l/4J u +1/4 J X2 + 3I2J 13 )mV 2J 13 
A 14 = l/4J u r + l/4J u s+(l/4J u +1/4 J 12 +3/2 J l3 )t + H2J X3 
A 1? =l/4J 12 r + (l/4J u +3/2 J l2 +l/4J u )s -1/2 J 12 +1/4 J l2 t 
A i m = (3/2 / u +1/4J 12 +l/4J u )r + l/2J u +1/4^ 35 + l/4/ u ? 
A xl 3 =l/4J lj2 r + (l/4J 11 + 3/2J l2 + l/4J l3 )s + l/2J l2 +l/4J l2 t 
A i16 = ( 3/2 J u + l/4/ u + l/4/ u )r-l/2/ u +l/4J 1 ^+l/4j 11 t 

A, 19 = (“ 4,3 “4,2 “3 1,1 ) r + (“ 4, 3 — ^ ^1,2 “ 4,1 ) s + (“ ^ -A, 3 “ A, 2 “ 

4.2 = 1 / 4 / 23/" + 1 / 4/ 23 5 + (l /4 J 21 + 1 / 4 / 2j2 + 3 1 2 J 2 ^t - M 2 J 2 3 

A 2 5=1/ 4/9 + 1 / 4 / 2 35 + (l/ 4 J 2 j + 1 / 4 / 2 2 + 3 / 2 J 2 3 )? + 1 / 2 / 2 3 

^9 g = 1/ 4/9 9/* + (l/ 4/9 j + 3/ 2/9 2 + 1 / 4/9 3 js — 1 / 2/ 2 2 + 1 / 4 / 2 9 ! 
/ 2 n = (3/2/2 i+l/ 4/ 2 2 + 1 / 4/9 3 j / ‘ + 1 / 2/ 2 j + 1 / 4 / 2 + 1 / 4 / 2 j / 

/ 2 14 = 1/ 4/9 2 /* + (l/ 4/9 1 + 3 / 2/9 2+1/ 4/9 3 + 1 / 2/ 2 2+1/ 4 J 2 9! 

4 , 17 = (3/2/2,! + 1 / 4 / 2 2 + 1 / 4 / 23 )r - 1 / 2/ 2 ! + 1 / 4 / 2>1 j + 1 / 4 / 2j / 

4.20 = (“ 4,3 “ 4,2 “ 3/ 2 i )/" + (- 4,3 ~ 3/ 2j 2 “ 4,1 ) 5 + (“ 3 / 2,3 “ 4,2 

4.3 = l/ 4 / 3>3 r + 1/4 /3 ! 35 +( 1/4 / 3i i +1/4/32 +3/2/33^-1/2/33 
4 >6 = 1/4/33^+1/4/335 +(1/4/3 ! +1/4/32 +3/2/33 )t +1/2/33 
4 9 = l/ 4 / 3j2 r + (l/ 4 / 31 + 3 / 2 J 32 +1/4/33)5-1/2/32 +1/4/32! 
4 >12 = (3/2/3! +1/4/32 +1/4/33 ) r + 1/2/34 +i/4/ 3>1 5 + 1/4/3/ 
4 i 5 = 1/4/32^ + (l/ 4 / 3jl +3/2/32 +1/4/33)5 + 1/2/32 +I/4/32! 
4 J8 = (3/2/31 +I/4/3 2 +I/4/3 3 )r -I/2/3 j +l/4/ 31 5 + I/4/3 / 

4.21 = (“ 4,3 “ 4,2 “ 3/3,1 )/" + (“ 4,3 “ 34,2 “ 4,1 ) 5 + (“ 34,3 “ 4,2 
/ 4j2 =1/4/3^ + 1/4/335 + (i/ 4 / 31 +1/4/32 +3/2/33)? -1/2/33 
4,3 =1/4/23/- + 1/4/2, 35 + (1/4/2,! + 1 / 44,2 +3/2/23)! -1/2/23 
4,5 = 1/4/33/- +1/4/335 +(1/4/3 ,! +1/4/33 +3/2/33)? +1/2/33 
ZI4 6 = 1 / 4 / 2 3/- + 1/4/335 + (l/ 4 / 2 i + 1 1 4 J 2 2 + 3/ 2/93)? + 1/ 2/33 
4 >8 =1/4/33/- + (1/4/33 +3/2/33 +1/4/33)5-1/2/33 +1/4/33? 
ZI4 9 =1/4/22/* + (1/4/91 + 3 / 2/ 2 2 + 1 4 / 2 3 ).s' — 1 2 J 2 2 + 1 / 4/ 2 2 ! 
/ 4 ii = (3/2/33 +1/4/33 + 1/4 4,3) r + 1/2 4 .i +1/4/335 + 1/4/33? 

ZI4 12 = (3 / 2 J 2 1 +1/ 4 / 2 2 +1/ 4 / 2 3)/* + l/2/ 2 1 +1/ 4/9 15 + 1/ 4 / 2 1? 

4,14 = 1 / 4 / 33 /- + ( 1 / 4/33 + 3 / 2/33 + 1 / 4 / 33 ) 5 + 1 / 2/33 + 1 / 4 / 33 ? 
ZI4 15 =1/4/99/- + (l/ 4 / 9 i + 3 / 2/92 +I/4/93 ) 5 +l/ 2 / 22 + I/4/22! 
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^ 4,17 

^4,18 

^4,20 

^4,21 

A,i 

^ 5,3 

^5,4 

A,e 

A,7 

A, 9 

A,10 

^5,12 

^5,13 

AlS 

A, 16 

^ 5,18 

^ 5,19 

^ 5,21 

Ai 

A, 2 

A, 4 

As 
Aj 
A,s 
A, io 
A, ii 
^6,13 
^ 6,14 
^ 6,16 
^ 6,17 
A,19 


H,20 


= (3/2 J 3l + \/4J 32 +\l4J 33 )r-\l2J xx +\l4J xx s + \l4J 3l t 
= (3/2J 2>1 +l/4J 22 + l/4J 23 )r-l/2J 2l +1/ 4J 2l s + 1/ 4J 2l t 
= ( — A , 3 ~A,2 “ 3^3,1 ) r + (“ ^3,3 ~^A,2 - A,1 ) 5 + ( - 3 / 33 ~ A,2 
= ( — A,3 ~A,2 ~^A,l \ + ( — J 2,3 _ 3/ 2 ,2 _ ^ 2 ,l ) 5 + ( _ 3/ 2 ,3 - ^ 2,2 
= \l4J 33 r + \l4J 33 s + {^l4J 3X + H 4J 32 +3 / 2J 33 )t -H 2J 33 
= l/4J l3 r + l/4J l3 s + {l/4J u +l/4J l2 +3/2J l3 )t-l/2J h3 
= l/4J X3 r + l/4J X3 s + (v4J 3A + \l4J X2 +3l2J 33 )t + \l2J 33 
= 1/ 4J l3 r + 1/ 4J l3 s + (l/ 4J xl + \l4J xl +3l2J l3 )t + \l2J l3 
= l/4J X2 r + (l/4J 3l +3/2J X2 +1/ 4J 33 )s -V2J 32 + 1/4 J 32 t 
= l/4/ u r + (l/4J u +3/2 J X2 + 1/4 Jj 3)5 -l/2J X2 + M4J X2 t 
= (3/2 J 3X +1/4 J 32 +l/4J X3 )r + l/2J 3l +V4J xl s+l/4J X it 
= (3/2 J xl +1/ 4 Jj 2 +l/4i 1>3 )r + l/2i u +l/4i 11 s + l/4J 1 ^ 

= l/4 A, 2 r + ( ll4 A,i + 3 / 2J 32 + 1 / 4 / 33)5 + 1 / 2 / 3 ^ + 1 / 4 / 32 ? 
= l/4/ 12 r + (l/4 / 1;1 +3/2 J X2 + 1 / 4 / 33)5 + 1/2 / 1j2 +1/4 J X2 t 
= ( 3 / 2/3 ^ + 1 / 4/3 2 + 1 / 4/3 3 ) 7 - -l/2 / 31 + 1 / 4/ 3j1 5 +I/ 4/3 / 

= (3/2/j j + 1/4/j 2 +l/4/ u )r-l/2 / 1-1 + 1 / 4 / u 5 + 1 / 4 / u ? 

= ( -< ^3,3 _ ^3,2 *^ 3 / 31 ) 7 - + (-/ 3; 3 “ 3/ 3 2 -^i)^ +(-3/ 3 _3 -/ 3 2 

= (“^1,3 “^1,2 - 3 / 1 j 1 )t- + (-/ 1j3 - 3 J X1 - J x j)s + (- 3/ 13 -/ 1>2 - 
= l / 4 J X3 r + l / 4 J X3 s + ( l / 4 J 2l + 1 / 4 / 2i2 + 3 / 2/2,3 )?-l/ 2/3,3 
= l/ 4 / u r + l/ 4/ 13 5 + (1/4/3,! +1/4/32 +3/2/33)? -1/2/33 
= 1/4/2,37- + 1/4/2,35 + (l/ 4 / 2 ,i + l/ 4 / 2 , 2 +3/2/23)? +I/2/2, 3 
= 1/4/3,37- + 1/4/3,35 + (1/4/33 + !/ 4 ^ i ,2 +3/2/33)? + 1/2/33 
= 1/4/2, 2 r + (i/ 4 / 2 , 1 + 3 / 2/ 2 ,2 +1/4/2, 3)5 -1/2/2, 2 +1/4/2, 2 1 
= l/ 4 /i 2 7 ‘ + (l/ 4 / u + 3 / 2 /i ,2 +1/4/3 3)5-1/2/3,2 +1/4/3 2? 

= ( 3 / 2/21 + 1 / 4 / 2, 2 + 1 / 4/ 2 ,3 )t" + 1 / 2 / 2 ,i + 1 / 4 / 2 , 35 + 1 / 4 / 2 , 3 ? 

= ( 3 / 2/11 + 1 / 4/12 + 1 / 4/33 )r + 1 / 2/33 + 1 / 4 / 3 , 35 + 1 / 4 / 33 ? 

= 1 / 4 / 2 2 7* + (l/ 4 / 2 3 + 3/ 2/-) 2 +1/ 4 / 2 3)5 + 1 / 2 /, 2 +1/ 4 J 2 2 t 
= l/ 4 /i, 2 t- + (1/4/33 + 3 / 2 /i, 2 +1/4/33)5+1/2/3,2 +1/4/32? 

= ( 3 / 2 / 2 ,! +1/4/2 , 2 +1/4/ 2 , 3 )t-- I/ 2 / 2, 1 + 1 / 4 / 2 , 35 + 1 / 4 / 2 , 3 ? 

= ( 3 / 2/33 + 1 / 4/32 + 1 / 4 / 33 ) 7 - - 1 / 2/33 +1/4/i,i5+1/4/i,i? 

= ( - «^2,3 — ^2,2 ^3/ 2 j )/- + (— *^2,3 ~3J 2 '2 ~ A,l) S + (“ 3^2,3 ~ A,2 
= ( -< A,3 ~J\,2 + ( -< A ,3 _ 3/i,2 +(-3^1,3 *^1,2 - 
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